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ABSTRACT 

The use of an adaptive mesh refinement 
procedure for analyzing high speed compressible flows 
with strong inviscid-viscous interactions is described. 
The adaptation procedure which uses both quadrilateral 
and triangular elements, is implemented with the multi- 
step Galerkin-Runge Kutta scheme. Elements that lie in 
regions of strong gradients are refined based on invisdd 
and viscous indicators to obtain better definition of flow 
features. The effectiveness of the finite element 
procedure is demonstrated by modelling Mach 11.7 flow 
over a 15 degree ramp. Numerical results are compared 
with predictions of the strong interaction theory and 
experimental data. The influence of factors such as 
mesh spacing and numerical dissipation on the 
accuracy of computed results is discussed. 

NOMENCLATURE 

A element area 

C Chapman-Rubesin parameter eq. (33) 

Cf skin friction coefficient 

Ch heat transfer coefficient, Stanton number 

Cp pressure coefficient 


c 1 ,...,c 4 Runge-Kutta constants 

d 1 • d 2 artifical dissipation fluxes 

Ei , E2 flux components 

e element error 


F a advective flux 


F v viscous flux 

heat flux 
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G growth factor, eq.(35) 

i unit vector defining velocity change, eq. (7) 

If . 12 components of unit normal vector 

M Mach number 

(MJ mass matrix 

[N] element interpolation function 

n unit vector normal to surface 

n x, n y components of unit vector n 

qn heat flux normal to surface 

p pressure 

Pr Prandtl number 

R global nodal residual 

r coordinate index , eq. (10) 

Re Reynolds number 

S solution domain boundary 

T temperature 

T x, Ty surface tractions 

t time 

t unit vector normal to i 

U typical conservation variable 

u, v velocity components in coordinate 

directions 

xi , X 2 coordinate directions 

Y maximum element error 

p key variable for refinement 

a threshold value for refinement 


1 



e element refinement indicator eq. (3) 

v Courant number 


p threshold value for derefinement, grid 

stretching parameter 

5jj Kronecker delta 

Y ratio of specific heats 

H2, P4 dissipation constants 

tw shearing stress 

X interaction parameter, eq. (33) 

Ax, Ay grid spacings in coordinate directions 


Subscripts 

aw adiabatic wail 

d dynamic conditions 

i summation index 

j index for nodes 

o total conditions 

res reservoir conditions 

w wall quantities 

oo freestream conditions 


Superscripts 

k index for multistep scheme 

n time step 


INTRODUCTION 

The potential role of computational fluid dynamics 
(CFD) in the design and analysis of high speed vehicles 
is being clearly defined by the national aerospace 
plane, which is envisaged to have applications at flight 
speeds exceeding Mach 25. Computational techniques 
that provide good understanding of the flow features and 
accurate estimates of aerothermal loads is essential in 
the design of such vehicles since ground based facilities 
cannot simulate the entire flight envelope. Description of 
the complex inviscid-viscous interactions with vehicle 
surfaces is possible with the full Navier-Stokes 
equations. At the Aerothermal Loads Branch at NASA 
Langley Research Center, finite element methodology 
is being developed for integrated fluid-thermal-structural 


analysis which can accurately predict the heating rates, 
pressure loads and the thermal/structural response. 

Compressible flows may contain discontinuities 
such as shocks as well as regions of high gradients 
such as boundary layers and shear layers which need to 
be adequately resolved. Since these flow phenomena 
occur over small distance scales which are not known 
apriori, the computational mesh needs to be adapted to 
model such high gradient regions without the use of an 
excessive number of elements in low gradient regions. 
Adaptive mesh refinement procedures can be used 
effectively to resolve details in the flow region and 
minimize elements elsewhere. Such procedures lead to 
highly unstructured meshes. Since finite element 
methods are characterized by their ease in handling 
completely unstructured meshes and their ability to 
include mesh refinement procedures these methods are 
pursued. Explicit finite element schemes have also 
demonstrated their capacity to produce good results for 
a variety of flow situations and configurations 1 ‘ 4 - 

The purpose of this paper is to describe the 
application of an adaptive refinement procedure to 
predict viscous compressible flow features. Adaptive 
mesh refinement procedures for compressible high 
speed flows are of recent origin. Mesh refinement 
procedures for triangular finite element meshes were 
initially detailed by Zienkewicz, Lohner, and Morgan^, 
and the application of these procedures to steady 6 , and 
transient' compressible flow problems has been 
demonstrated extensively. Adaptive procedures for 
finite element meshes with quadrilateral elements have 
been developed by Oden et. al 8 and by Shapiro and 
Murman 8 A mesh regeneration scheme developed by 
Peraire, et. al. 10 has found application in generating 
meshes comprised of triangles (in 2D) and tetrahedrons 
(in 3D) for inviscid flows. 

Most of the research in adaptive finite element 
methods has been in inviscid flows. Application of 
adaptive procedures to model complicated high speed 
viscous flows are virtually non-existent. A notable 
exception has been the recent effort to use an upwind 
finite element formulation 1 1 and unstructured triangular 
grids to detail complex flows resulting from shock 
interaction on a cylindrical leading edge. 

The present study extends the adaptive 
refinement procedure and finite element formulation of 
reference 12 to predict inviscid and viscous flow details 
for hypersonic flows with strong inviscid-viscous 
interactions. 

MESH ADAPTATION 

The classical finite element mesh refinement 
scheme is the addition of elements in regions of high 
gradients. Elements that lie in these regions are divided 
into smaller ones by a subdivision process. Both 
triangular and quadrilateral elements can be enriched 
by adding a central node and/or midside nodes. Figure 
1(a) illustrates a triangular mesh that has been locally 
enriched to capture an oblique shock. For a typical 
quadrilateral element such a subdivision can result in 
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( 4 ) 


the generation of four smaller elements with the 
possibility of the presence of midside nodes not being 
connected to neighboring elements. These midside 
nodes are sometimes called ’hanging nodes.’ Figure 
1 (b) illustrates a mesh with hanging nodes. The solution 
algorithm must be modified to account for hanging 
nodes, typically by introducing constraint equations. 
One way to avoid these unconnected nodes is to 
transition from a crude quad mesh to a fine one using 
triangular elements. Figure 1(c) shows this type of 
refinement which needs no special constraint equations 
and the procedure provides a ’natural" way to transition 
between meshes of different density. Automation of the 
adaptive procedure is accomplished with refinement 
indicators. 

Refinement Indicators 

The rationale for using refinement indicators is 
that, while it is possible to predict the location and 
strength of shocks, boundary layers, etc. for some simple 
flow situations, the analyst in general will not have prior 
knowledge of the location of regions containing sharp 
changes in flow variables. The decision to refine a 
particular region of the mesh can be based on either 
a-priori or a-posteriori error estimates. The procedure 
adopted in this paper is to complete an analysis on a 
given mesh and then refine the mesh at a certain stage 
in the analysis. The refinement at this stage is based on 
inviscid and viscous "error" indicators computed on the 
initial mesh. The new mesh is used for the subsequent 
analysis and then further refinements are performed if 
needed. 

The aim of the adaptive mesh refinement is the 
minimization of maximum errors that occur in the finite 
element domain. The mesh requirement is then 

minimize (max ej ) 1 <i<n (1) 

where i denotes the element number, n the number of 
elements in the domain and ej the error in the i-th 
element. An optimal mesh then satisfies the equal 
distribution condition, 

e. = constant (2) 

which indicates that the error measure is uniformly 
distributed for all elements of the finite element mesh. 
The adaptive mesh refinement procedure is designed to 
add elements in appropriate regions to satisfy Eq. (2). 
Since the solution is not known a-priori, an 
approximation to the error measure is computed using 
the finite element solution. If 'p* is considered the key 
variable representative of overall solution behavior, then 
an error measure for invisdd refinement is given by 

j.k- 1.2 

where 0jk are scaled second derivatives with a typical 
term 612 defined as 


0 i 2 " e *y " P-*y AUlpJ + ej^ lp.,1,^1] 

ui 

Here nd is the number of nodes in element i, and the 
comma implies partial differentiation. The first term in the 
denominator of equation (4) scales the second 
derivative of the key variable to ensure that shocks of 
different strengths are adapted equally. This prevents 
the strongest shock from being overdefined at the 
expense of capturing weaker shocks. The coefficient e in 
the second term is used to smooth the indicator in 
regions of oscillations. 

The procedure adopted in this paper is to refine 
all elements that satisfy the criterion 

e, > aY (5) 

and derefine all elements that satisfy 

e,<pY (6) 

where a and p are preset threshold constants and Y the 
maximum element error over the entire domain. The key 
variable used for the inviscid refinement is typically 
pressure or density and a and p are usually 0.2 and 0.1 
respectively. 

The refinement indicator for viscous regions is 
based on a variation of the explicit dissipation scheme 
detailed in reference 13. It uses the notion that near the 
wall the velocity changes within the boundary layer are 
maximum in the direction normal to the predominant 
velocity component. Unit vector i defines the change in 
the value of the velocity vector V and is given by, 

ij-yyv.f + v.? M 

where V 2 = (u 2 + v 2 ). The expression for element error 
is given by 

e ( = h 2 |^ (V.t)l (8) 

where unit vector t is defined normal to i and is 
illustrated in figure 2. As with the inviscid indicator, 
elements with errors above a preset tolerance value are 
refined while those with errors below the threshold are 
derefined. 

Adaptive Procedure 

The starting point of the adaptive procedure is a 
’skeleton’ mesh or base mesh which contains only a few 
elements. The skeleton mesh, which consists of 
quadrilateral elements, needs to be completely 
structured to begin the adaptive procedure. A structured 
mesh implies that each interior node in the domain is 
surrounded by the same number of elements. Prior to 
the first analysis, the mesh refinement program does an 
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overall refinement where each element of the ’skeleton* 
mesh is subdivided into four quad elements. This 
procedure is repeated a few times until the nodal density 
of the resulting mesh is deemed sufficient to obtain a 
solution which captures the main details of the flow field. 

Refinement indicators are computed based on the 
solution obtained on this ’initial* mesh and elements that 
need to be refined or derefined are identified. All 
elements in the mesh that have indicators above the 
preset refinement threshold value are enriched while 
those elements that have values below the threshold 
derefinement value are coarsened. The refinement 
strategy used is such that at each mesh change, only 
one level of refinement or derefinement is permitted. On 
refinement of a typical element, the "sub-elements" that 
result could be all quads, or a combination of quads and 
triangles. The number and type of the resulting ’sub- 
elements’ depends on the refinement level of elements 
that surround this element. 

Figure 3 shows the elements that result in a 
typical refinement and coarsening procedure. Figure 
3(a) shows the original mesh where elements B, C, and 
D are to be refined. The mesh that results is seen in Fig. 
3(b). If on this refined mesh, element group C, which 
includes "sub-elements’ Cl , C2, C3, and C4, needs to 
be coarsened, the mesh that results after derefinement 
appears in Fig. 3(c). 

Boundary Layer Refinement 

In addition to the refinement based on error 
indicators detailed in earlier sections, the facility to 
generate layers of structured quad elements at the wall 
is also included within the framework of the analysis 
program. This technique is developed to ensure a high 
density of nodes at the no-slip boundaries to resolve 
the thin laminar hypersonic boundary layer. The 
program divides the elements located at the wall into 
a specified number of sub-elements and elements thus 
generated can be refined and coarsened at successive 
refinement levels. The nodal coordinates of these new 
elements are defined such that the mesh generated is 
stretched normal to the wall. The location of the nodes 
along normal lines at each wall location is obtained 
from the expression of Roberts 14 , 


. (P+1) - 0-1) / [{(P+t) / 0-i)} lH l , „ 

O. = I < p < oo 

1 I(P+1)/(P-1)1 +1 

x=x* + Ax8 (9) 

j = 1 , nnew 

y=y w + Ay5i 

where nnew is the number of new elements that result 
from the subdivision. Ax and Ay the grid spacings in the 
coordinate direction, p the grid stretching parameter, 
and r the coordinate index given by, 

r = j/nnew (10) 


where j is the index locating the node from the wall in 
the normal direction. The nodal unknowns for the new 
grid points at each wall location is linearly interpolated 
from the nodal values of the parent element. The 
meshing procedure clusters more nodes near the wall 
as p approaches unity. 


FINITE ELEMENT FORMULATION 

The explicit multistep Galerkin-Runge Kutta 
formulation used in reference 12 for inviscid flow is 
extended to model the compressible Navier-Stokes 
equations. The solution algorithm is applied to the flow 
equations in conservation form 


U*t + Ej i = 0 (11) 

where U is the vector of variables and Ej the flux 
vectors. The flux vectors can be written as 


E, = F a -F"-F h d2) 

i i i i 

where F a , F 7 and F h represent the advective, viscous 
and heat fluxes respectively. The vectors of conservation 
variables and fluxes are given by 



where p is the density, Uj the velocity components in the 
coordinate directions, p the thermodynamic pressure, e 
the total energy, qj the heat fluxes, xjj the viscous stress 
components and 8y the Kronecker delta. The shear 
stresses and heat fluxes are given by, 


vKkV^i+V 

q-'kT, 


(14) 


where p and X are the coefficients of viscosity, k the 
thermal conductivity and T the temperature. Other 
constitutive relations employed include: 

e = C v T + (i i 2 +v 2 )/2 (15) 

p = pRT 


4 



where Cv is the specific heat at constant volume and R 
the gas constant. The coefficients of viscosity are related 
by Stokes hypothesis 


and D2 and O4 are second and fourth difference 
dissipation terms which may be 'frozen* for each 
timestep. 


X = "§^ (16) 

and p is assumed a function of temperature and 
obtained by Sutherland's law. The thermal conductivity 
is computed from the Prandtl number Pr which is 
assumed constant along with constant specific heats. 
Equation (11) is solved subject to proper initial and 
boundary conditions which include: (a) specification of 
conservation variables on the inflow plane, (b) 
imposition of no-slip conditions and specified 
temperatures on aerodynamic surfaces, and (c) outflow 
surface integrals provided by the finite element 
formulation. 


Galerkin-Runge Kutta Algorithm 

The system of conservation laws shown in Eq. 
(11) can be written as, 

U t = -E|,| (17) 

Application of the method of weighted residuals results 
in the finite element equations, 

J[N]{^dA = -J[N] (E. .)dA (18) 

A A 

where [N] is the element weighting and interpolation 
function and A the domain of interest. Integrating by 
parts on the RHS terms yields, 


The second difference artificial dissipation term is 
needed to stabilize solutions in the presence of sharp 
gradients. This dissipation operator is given by, 


and 

d. = KU,. ( i not summed ) (22) 

where kj control the amount of dissipation added and 
are given by, 

p.h 4 

•S = -fr- I P.jj I ( i not summed ) (23) 

where p2 is the second difference dissipation constant 
and p the pressure. The first derivatives of pressure are 
calculated as element quantities and nodal second 
derivatives are then obtained from the relation, 

jNp >xx dA = -Jp^N, x dA (24) 

A A 

where the comma implies differentation. 

The linear fourth difference operator is needed to 
damp out spurious oscillations and provide background 
dissipation and is given by, 

D 4 = ^4 (AX i )4 U -iiii (25) 


[M] {^-} = J{N,.} [N] dA {E} -J{N) [N] dS [l.E J (19) 

A S 


where the fluxes Ej are interpolated the same way as the 
conservation variables. 


The fourth differences are evaluated by repeated 
application of a 9-point Laplacian stencil. Typical values 
of the dissipation constants P2 and P4 are 1/10 and 
1/200 respectively. 

Calculation of Wall Coefficients 


The time marching scheme is similar to the multi- 
step time integration scheme which appears in 
reference 9 and can be written as, 


U <k) =U n + c k [^- R(U M ) + D]k=1,..,4 
jj 

u n+1 =u (4) 

i i 


( 20 ) 


where the dissipation operator D can be written as, 


D 


_n 

D + D 

2 4 


( 21 ) 


Heat transfer, skin friction and pressure 
coefficients are essential for design considerations since 
they provide a direct measure of the severity of 
aerodynamic loads on flight surfaces. Experimental 
measurements for hypersonic flow behavior yield 
surface distributions for pressure loads, heat transfer 
rates and skin friction coefficients. Accurate computation 
of the wall coefficients in the finite element 
methodology is essential since their comparison with 
experimental data forms the basis for code validation. 

The heat transfer and shear stress coefficients 
are related to the gradients of temperature and fluid 
velocity at the surface. With a structured computational 
grid it is possible to calculate derivatives at the wall by 
simple differencing techniques and obtain results of 
desired accuracy. For completely unstructured finite 
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element meshes the calculation of wall derivatives by 
differencing is rather complicated. Computation of heat 
fluxes and shear stresses at the wall is possible in the 
finite element context by using the ’consistent’ 
calculations. These calculations, introduced in 
reference 13, use equation (19) as the basis for 
computing wall coefficients. 

The finite element equations (eq. 19) can be written as 

M - f< N > M dS l l i E J ( 26 ) 

A S 


At steady state the transient terms approach zero, 
implying 

J{N) [N] dS [I.EJ =J{N, j ) IN] dA (E) (27) 

S A 

The x and y momentum equations for elements on the 
surface reduces to. 


C H - q«/[P-U- I(e + p/ p + u%)- - (e + p/ p )J 

C, = VP~ U ~ (32) 

C p = Pv/(jP- u ^) 

where p„ , p„ and u„ are freestream values and p* is 
the wall pressure. 

COMPUTATIONAL RESULTS 

The ability of the adaptive finite element 
procedure to predict flow details for viscous flows is 
illustrated by modelling the hypersonic flow over a ramp. 
This problem has special interest due to its application 
in the design of control surfaces for high speed vehicles 
such as the aerospace plane. Theoretical estimates and 
experimental results are available to gauge the 
accuracy of the finite element calculations. 



[N] dS {T x } = (R x ) 


1 


(N) [N]dS{T y } = {R y ) 


(28) 


where T x and Ty are the surface tractions and R x and 
Ry the residuals for an element. These equations are 
assembled for all the elements that lie on the surface. 
Diagonalizing the matrices defined by the surface 
integrals in equation (28) yields an explicit procedure for 
obtaining nodal values for the surface tractions. The 
surface tractions are related to the stress tensor by the 
relation, 



where xq are components of the stress tensor and n x 
and n y are components of the unit vector n, normal to 
the surface. The wall shearing stress is then given by, 


x w = T y n x -T x ny (30) 

Considering the energy equation for elements that lie on 
the surface, at steady state the finite element equations 
can be written as, 

J{N) [N] dS { qn } = (FW (31) 

s 


where q n is the nodal heat flux for an element. 
Diagonalizing the coefficient matrix and assembling the 
element equations yields nodal values of the heat flux 
on the surface. The coefficients of heat transfer (Ch). 
skin friction (Cf) , and pressure (Cp) are computed from 
the relations. 


Flow Description 

The problem considered is illustrated 
schematically in Figure 4. Inflow at Mach 1 1 .68 interacts 
with the leading edge producing a shock due to the 
boundary layer displacement effect. The boundary layer 
does not have enough momentum to overcome the 
adverse pressure gradient generated by the ramp 
resulting in flow separation and recirculation near the 
corner. The separated boundary layer then reattaches 
downstream of the corner with the surface pressure 
rising through the separated and reattachment regions. 
The compression fan generated in the separation zone 
eventually coalesces downstream to form the induced 
shock. This shock interacts with the leading edge shock 
to produce a resultant shock, an expansion fan and a 
shear layer or slip line. 

The interaction between the inner viscous shear 
layer and the outer inviscid flow can be categorized 
based on the nature of coupling. At the leading edge of 
the plate the induced shock and the boundary layer are 
essentially of the same order and the effects of the 
interaction can be reasonably predicted by various 
interaction theories. At the corner the separation results 
upstream of the induced shock due to the coupling 
between the shock and the thickening of the boundary 
layer due to compression. This type of interaction 
necessitates the need for solution of the full Navier- 
Stokes equations. 

Strong Interaction Theory 

Near the sharp leading edge the displacement 
thickness increases from zero and the flow turning 
outside of the boundary layer causes compression and 
a leading-edge shock. The strength of the shock 
depends on the incident Mach number and has a strong 
effect on the growth of the boundary layer at the leading 
edge. Very near the leading edge, the pressure at the 
edge of the boundary layer is very high compared to the 
freestream pressure marking the region of ’strong* 
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Experimental Measurements 


interaction. Downstream of the leading edge the 
interaction between the inviscid and viscous flow 
’weakens' and the wall pressure expands to freestream. 

From laminar boundary layer theory the 
appropriate similarity parameter for the shock-boundary 
layer interaction is the Lees-Probstein parameter X 
given by 

X = ut/C/Jte x (33) 

where C is the Chapman-Rubesin compressibility 
parameter, M„ the freestream Mach number, and Rex 
the local Reynolds number. The effects of the coupling 
between the leading-edge shock and the boundary 
layer have been predicted by various strong 15,16 and 
weakl 7 interaction theories. The complete interaction 
theory of Bertram and Blackstockl 3 unifies the strong 
and weak interaction theories by bridging their limits of 
applicability. The assumptions used in the development 
of these interaction theories include: 

(a) the flow is hypersonic (Moo »1) 

(b) no gas dissociation 

(c) Prandtl number constant 

(d) vorticity effects due to the curved leading edge 
shock are negligible, and 

(e) pressure distibution is assumed to vary 
exponentially with distance from the leading 
edge. 

The complete theory relates the surface pressure 
to the interaction parameter X as follows, 

= 0.83 +|- GX (34) 

where G is the growth factor defined as, 

G = 1 .648 + 0.352] (35) 

2 aw 

and Taw is the adiabatic wall temperature given by 

T aw = T »M n£] (36) 

for laminar flows. The variation of heat transfer 
coefficient with surface temperature in the strong- 
interaction region is tabulated by Li and Nagamatsu 1 ® . 
For y = 1 .4, variation of Ch with wall temperature T w is 
given by 

C H [Re* /X] 1/2 = C 1 ( 37 ) 

where the constant Ci depends on T w and the ratio of 
specific heats and is obtained from reference 19. 


The experimental studies used for validation of 
the finite element approach were conducted in the 
Calspan 48-inch shock tunnel. These experiments 20 
have generated heat transfer, surface pressure and skin 
friction measurements for hypersonic conditions ranging 
from Mach 11.5 to Mach 1 8.9. 

The Calspan shock tunnel consists of a shock 
tube driver, followed by a conical expansion section and 
a contoured nozzle to produce uniform axial flow in the 
test section. The tunnel is started by rupturing a double 
diaphragm which causes expansion of the air in the high 
pressure driver section. A normal shock is generated 
with the high temperature and high pressure air 
between the normal shock and the driver-driven 
interface. When the shock front strikes the end of the 
driven section it leaves a region of stagnant high 
pressure air. This air is then expanded through a 
contoured nozzle to desired freestream conditions. 

Skin friction transducers, thin-film resistance 
thermometers and pressure transducers mounted at 
suitable locations obtain skin friction, heat transfer, and 
surface pressure measurements. Details of the 
experimental setup are found in reference 20. 

The flow parameters used in the numerical 
computations correspond to the nominal test conditions 
for hypersonic flow over a ramp presented in reference 
20 and given below 

Moo =11.68 

Re«o = 1 .71 4E+05/ft 

Tres - 2989 R 

Pd = 0.3589 psia 

The experimental data to be used for comparisons will 
consist of measurements obtained for two sets of tunnel 
runs^l . These measurements correspond to runs 19 
and 21 , and flow conditions for these runs are as follows 

Run 19 Run 21 


Moo 

= 11.7 

Moo 

= 11.64 

Reoo 

= 1 .734E+05/ft 

Reoo 

= 1.644E+05/ft 

Tres 

= 2911 R 

Tres 

= 3105 R 

Pd 

= 0.3458 psia 

Pd 

= 0.369 psia 


The measurements of runs 19 and 21 yield surface 
heating rates, pressures and skin-friction. The 
coefficients of pressure and skin-friction are obtained by 
dividing experimental values by the freestream dynamic 
pressure. The coefficient of heat-transfer is obtained 
from the heat flux at the wall by the relation, 

C H = P.u. Y cXs- t J (38) 

where T w , the wall temperature, assumed to be the 
same for both runs, is obtained from reference 20 and 
the specific heat is assumed constant. 
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Computational Strategy 

The hypersonic flow over a ramp results in 
invisdd-viscous interactions that have distinctly different 
character in various regions of the domain. Near the 
leading edge, the mesh needs to be well refined and 
uniform to capture the gradients in both the streamwise 
and normal directions. Away from the leading edge, the 
flow is basically the development of the viscous 
boundary layer and mesh spacings at the surface in the 
normal direction are of concern. At the corner, the 
coupling between the inviscid and viscous regions 
necessitate the need for adequate refinement to capture 
the separated and attached regions of flow. 

The adaptive finite element approach was 
implemented on the NAS Cray-2 supercomputer. The 
analysis is started on an initial mesh and iterations are 
performed until the L2 norm of the residuals of the 
conservation variables dropped over three orders of 
magnitude. The mesh is then refined/coarsened based 
on the inviscid and viscous error indicators. On the final 
mesh, convergence is monitored by observing changes 
in the residuals of the conservation variables, the heat 
flux on the surface and variations in pressure values 
near the surface at the outflow plane. A local time- 
stepping procedure that is based on the CFL and 
viscous stability requirements is used to accelerate 
convergence to steady-state and details of this 
procedure are found in reference 13. Typical 
computational speeds on the scalar code were about 
2.E-04 CPUs/timestep/node and vectorization efforts 
are underway to enhance this speed. 

An initial analysis of the complete flow domain 
resulted in inaccurate descriptions of the heating rates 
and surface pressures. This analysis indicated the need 
for the finite element mesh to have more elements at 
critical locations of the flow, such as, at the sharp 
leading edge and the corner to obtain satisfactory 
correlation with experimental data. The analysis 
described herein divides the flow domain into three 
regions. The three regions are shown in figure 5 and 
are categorized as : (a) the sharp leading edge section, 
(b) the flat plate section, and, (c) the ramp section. 

Sharp Leading Edge 

The finite element mesh that results after multiple 
refinements/derefinements depends a greal deal on the 
'skeleton* mesh used to start the adaptation procedure. 
A judicious selection of the skeleton mesh models the 
physics of the flow better and minimizes the number of 
elements used in the analysis. The skeleton mesh used 
to model the sharp leading edge is shown in Fig. 6. The 
skeleton mesh was designed such that on refinement, 
aspect ratios of unity exist right at the leading edge and 
higher aspect ratios prevail downstream of the leading 
edge. The mesh obtained after three refinements on the 
initial mesh is shown in Fig. 7(a) and contains 7245 
nodes and 7235 elements (5781 quads and 1454 
triangles). Density contours for the leading edge section 
of mesh in figure 7(b) show good definition of the 
leading edge shock and smooth variation of density 
within the boundary layer. At the leading edge, the 
boundary layer displacement effect causes the density 


at the wall to increase by an order of magnitude and 
then to drop off quite rapidly away from the leading 
edge. 

Details of trie finite element mesh very near the 
sharp leading edge ( box A in figure 7(a)) is shown in 
figure 8(a) and the density contours on this section of the 
mesh appear in figure 8(b). It is interesting to note the 
'wiggles* that appear in the contours at x=0.005 where 
the element aspect ratio jumps by a factor of 20 in the 
flow direction. The effect of this jump is also seen in 
figure 9 in the distribution of the wall pressures near the 
leading edge. 

A region of the mesh where the element aspect 
ratios change considerably normal to the flow direction 
is at y = 0.005. This jump is inherited from the macro- 
elements that constitute the skeleton mesh in figure 6. 
The effect of the aspect ratio jump in the normal direction 
is best seen by focusing in on the region at the exit 
plane (x=0.1) defined by box B in figure 7(a). The finite 
element mesh at the outflow plane is shown in figure 
10(a) and the distribution of pressure along y at x=0.1 is 
plotted in figure 10(b). The presence of oscillations near 
the wall and at mesh transitions is very evident in this 
distribution. The bigger 'wiggle* occurs where the 
element aspect ratio in the normal direction changes 
abruptly by a factor of ten. The analysis was repeated 
with both the second and fourth order dissipation turned 
off for elements located close to the wall and away from 
the vicinity of the sharp leading edge. This was possible 
using arrays that contain information about element 
neighbors. The distribution of pressure at the exit with no 
explicit dissipation close to the wall in shown in figure 
10(c). The pressure jump at the wall is absent and 
pressure levels within the boundary layer are higher. 
However, the oscillations caused by the aspect ratio 
jumps are still present. The distribution of pressure at the 
exit, shown in figure 11, indicates the presence of 
pressure gradients normal to the wall even at regions 
close to the wall. Pressure, being derived from the 
conservation variables, appears to be more sensitive to 
mesh transitions than, say, density, the distribution of 
which is plotted in figure 12, exhibiting a smaller kink at 
y-0.005. 

The surface pressure and heat transfer coefficient 
obtained from this analysis are compared to predictions 
by the interaction theories in figure 13. The results 
compare favorably, especially the heat transfer 
coefficient. Pressure levels close to the leading edge 
are well behaved and do not exhibit the 'inverse spike* 
reported by other researchers22. The maximum 
deviation in pressure levels is seen to be near the 
leading edge. 

It is instructive to take a quick look at the flow 
physics close to the leading edge. As seen in Figure 14, 
very close to the leading edge the flow is near free 
molecular flow and can be modelled only on the basis of 
kinetic theory. The transitional layer has mixed kinetic 
and continuum properties. The merged layer which lies 
a little downstream of the transition layer can be 
described by the Navier-Stokes equations but needs 
velocity slip and temperature jump conditions at the 
wal|23. Downstream of the merged layer, where a 
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distinct inviscid layer develops between the shock and 
the viscous layer, is the interaction region. The Navier- 
Stokes equations with no-slip boundary conditions do 
not accurately model the physics of the flow up to this 
interaction region. So the accuracy of the finite element 
calculations very close to the leading edge is less than 
perfect and discrepancies on comparison with the 
strong interaction theories is not unexpected. 

Flat Plate Section 

Profiles of the conservation variables obtained at 
the exit of the leading edge section are used as inflow 
conditions to the flat plate section. 

The skeleton mesh used to model the second 
section of the problem, the flat plate region, is shown in 
figure 15. After three overall refinements the elements at 
the wall were divided into 10 sub-elements according to 
equation (9). The grid parameter 3 was set to be 1.02 
and the subdivision yields elements with aspect ratios 
over 200 at the wall. The finite element mesh that results 
after three refinements on this initial mesh contains 8483 
nodes and 8593 elements (5151 quads and 3442 
triangles) and is shown in figure 16(a). Density contours 
on this mesh in figure 16(b) show the mesh capturing 
the leading edge shock and growth of the boundary 
layer along the plate. Details of the finite element mesh 
near the inflow (box A in figure 16(a)) and the density 
contours on this section of the mesh appear in figures 
17(a) and 17(b). Density profiles at the inflow and the 
exit of the flat plate section, shown in figure 18, illustrate 
growth of the boundary layer and weakening of the 
leading edge shock as the flow traverses the plate. 

A comparison of the density profiles at the 
outflow of the leading edge and at inflow to the flat plate 
section in figure 19 indicates the adequacy of nodal 
distribution in the normal direction. Linear interpolation 
is used at the inflow of the flat plate section which tends 
to dampen some of the oscillations present at the exit of 
the leading edge section as illustrated by the pressure 
profiles in figure 20. It is to be noted that these 
oscillations become fixed at the inflow planes, and 
hence will influence subsequent analyses. 

As with the leading edge section, profiles on the 
exit plane (x=0.7) of the flat plate section are used to 
define the inflow boundary of the ramp section. 

Ramp Section 

The skeleton mesh used for modelling the ramp 
section of the flow domain is shown in figure 21. After 
three overall refinements on this skeleton mesh a layer 
of structured quadrilateral elements were generated at 
the wall. Each element at the wall was divided into ten 
sub-elements with a stretching factor of 1.03. The mesh 
that results after three refinements on this initial mesh is 
shown in figure 22(a). Density contours obtained on this 
mesh appear in figure 22(b) and salient characteristics 
of the flow such as the leading edge, induced and 
resultant shocks are all seen to be well captured. The 
ability of the adaptive refinement procedure to locate 
regions of high gradients is exhibited by taking a closer 
look at the mesh in the ramp comer (box A in figure 


22(a)). The mesh in figure 23(a) and density contours 
on this mesh in figure 23(b) illustrate capture of the 
leading edge shock and compression of flow that results 
in the induced shock downstream of this region. The fine 
mesh at the surface is needed to model separation and 
subsequent reattachment of the boundary layer in the 
vicinity of the corner. Explicit dissipation terms for 
elements that are located close to the surface were 
zeroed out to ensure accurate computation of surface 
quantities. 

Details of the flowfield, especially the coupling 
between the inviscid and viscous regions of flow, are 
best described by profiles of select variables at specified 
axial locations. Figure 24 indicates the various x-stations 
wherein the variation in y direction of quantities such as 
pressure and temperature are studied. Locations 
indicated by A, B, C and D are used to describe flow 
details in the domain up to the comer while locations E 
to H detail flow characteristics up the ramp. 

The distribution of pressure normal to the surface 
at locations A to D is shown in figure 25. As the flow 
moves from the inlet toward the corner, the presence of 
the adverse pressure gradient due to the ramp is clearly 
seen. Fluid in the vicinity of the comer cannot overcome 
the effects of the pressure gradient as well as the skin 
friction at the wall and flow separation results. The 
separated boundary layer becomes a free shear layer 
outside of the recirculation region. Pressure profiles in 
figure 25 and velocity profiles in figure 26 illustrate this 
effect. The boundary layer separation is clearly seen by 
reversal in flow direction of the u-velocity profiles. The 
distribution of temperature at various x-stations in figure 
27 shows evidence of increase in boundary layer 
thickness as the flow approaches the ramp corner. 
Gradients of the temperature profiles at the surface 
indicate the decrease in heat transfer rates from location 
A to location D. 

Characteristics of the flow up the ramp are 
described by the distributions at locations E to H. 
Pressure profiles at E, F, and G in figure 28 indicate the 
strengthening of the induced shock and the change in 
location of the leading edge shock relative to the 
surface. The distribution of pressure at exit of the ramp 
(location H) shows the resultant shock which is derived 
from the interaction of the leading edge shock with the 
induced shock. A weak expansion fan also results from 
this interaction, the effect of which is seen by the small 
drop in surface pressures at location H relative to that 
at G. The distribution of density in figure 29 indicates the 
thinning of the boundary layer as the flow moves up the 
ramp. This reduction in thickness of the boundary layer 
is also seen in figure 30 which shows u-velocity profiles 
at various x-stations. The presence of the recirculation 
is seen by the flow reversal close to the surface at 
station E. At the exit from the ramp section, the boundary 
layer profile is fully developed and downstream of this 
region the interaction between the inviscid and viscous 
regions of the flowfield is negligible. Temperature 
profiles at locations up the ramp, shown in figure 31 , 
reflect the "thinning" of the boundary layer. The largest 
gradient in temperature at the surface is seen to be at 
location G indicating the approximate location of peak 
heating rates. 
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Heat flux, skin friction and pressure coefficients 
for the ramp surface obtained from the finite element 
approach are compared to experimental observations 
of Holden^ 1 . As mentioned earlier, flow conditions used 
in the finite element analysis were conditions 
representative of runs 19 and 21 in the Calspan tunnel. 
Figure 32 compares the finite element predictions for 
heat transfer coefficient to those obtained from tunnel 
runs. Heat transfer coefficients computed from the finite 
element scheme compare well with both the 
experimental runs for locations upto the corner as well 
as halfway up the ramp. Heating rates measured near 
the exit of the ramp for run 19 are about twice that of 
both run 21 and numerical predictions. The finite 
element predictions seem to compare better with 
experimental data of run 21 as seen by figure 33. The 
heat transfer coefficients match quite well at the 
separation and reattachment regions and the values of 
the peak heating rate are in excellent agreement. The 
main discrepancies seem to be location of the peak heat 
rate, which for the analysis is shifted to the left, and a 
steeper drop in heating rates downstream of this 
location. 

Skin friction coefficients calculated from the finite 
element approach are compared to those measured in 
tunnel runs 19 and 21 in figure 34. Levels of the skin 
friction coefficient compare well with the experimental 
data for regions unaffected by the recirculation while the 
size of the recirculation region is underpredicted by the 
finite element approach. Both experimental runs 
measure flow separation at around x=15" while the finite 
element predictions locate separation at x=16". The 
location of the reattachment point from run 21 is x=20.3* 
and the finite element procedure locates this point at 
x=19.3". Unsteadiness of the recirculation bubble, which 
could account for the discrepancies in the finite element 
predictions, was not considered in this analysis. 

A comparison of surface pressure coefficients 
using the finite element approach with data from 
experimental runs 19 and 21 appears in figure 35. 
Pressure levels predicted by the analysis is seen to be 
lower at the inflow and recirculation regions. The rise in 
pressure up the ramp is much more rapid with the final 
surface pressure falling midway between pressure 
levels of runs 19 and 21. The experimental data shows 
sharp peaks in pressure coefficients near the exit while 
the finite element approach indicates a gradual 
downward trend. This implies that the interaction 
between the leading edge shock and induced shock is 
stronger than that predicted by the finite element 
approach. The pressure rise for Mach 11.7 inviscid flow 
over a 15° ramp is about 17.5 and the finite element 
analysis yields a peak pressure rise of over 22. 

Refinement Indicators Revisited 

The finite element meshes for analyzing flow over 
the ramp were obtained by multiple refinements on the 
initial meshes, using the inviscid indicator, with pressure 
as the key variable, and boundary layer indicator. 
Proper selections of the key variable and values of a 
and 3, the refinement and derefinement thresholds, 


impact heavily on the quality of the resulting finite 
element mesh. The effect of refinement and threshold 
parameters on the finite element mesh for the ramp 
section merits investigation. Meshes that result from 
refinements are evaluated by comparing nodal 
distribution in the y-direction for some specified x- 
station. 

Four cases of interest, obtained by varying the 
adaptation parameters are described below. 

Case 1 : Key variable - pressure, a=.15, p=.l , viscous 
refinement indicator off. 

Case 2: Key variable - pressure, a=.15, p=.1, viscous 
a=.15. 

Case 3: Key variable - density, a=.15, p=.1, viscous 
refinement indicator off. 

Case 4: Key variable - density, a,p=,1, viscous 
refinement indicator off. 


The pressure profile at location A (figure 24) in 
the flat portion of the ramp for the initial mesh is seen in 
figure 36. For both case 1 and case 2, the use of 
pressure as the key variable causes the leading edge 
shock to be refined. The profile reflecting case 2 also 
shows nodes being added inside the boundary layer. 
The location of these new nodes are seen in figure 37 
from the u-velocity profiles. The viscous indicator is seen 
to add elements where the velocity profile is essentially 
linear. 

The use of density as the key variable (case 3) 
indicated by the profile in figure 36 shows refinements at 
the shock and the boundary layer. Elements within the 
boundary layer are refined due to sharp changes in the 
density profile as seen in figure 38. The effect of 
lowering the refinement threshold is shown in figure 39 
where the pressure profile that results from case 4 is 
compared to that of case 3. More elements at the 
leading edge shock are refined for case 4 and the 
pressure profile is very similar to that obtained using 
pressure as the key variable and the viscous indicator 
(case 2) in figure 36. 

The profiles in figures 36-39 illustrate the 
inadequacy of using equation (4) with pressure as the 
key variable to capture details of the boundary layer. 
The use of density as the key variable with appropriate 
threshold values causes refinements at the shock as 
well as within the boundary layer. The use of the viscous 
indicator equation (8) leads to excessive refinements 
within the boundary layer, where the velocity profile is 
essentially linear. 

Oscillations in pressure close to the wall is 
observed in figures 36 and 39. The pressure levels of 
the added nodes are higher than the values at the 
adjacent nodes. The inaccuracy is due to pressure 
being derived from conservation variables while 
conservation variables at a new node are obtained by 
averaging variables at adjacent nodes. A time history of 
pressure at such a typical node in figure 40 shows the 
nodal pressure dropping down to the right levels within 
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200 iterations on this mesh. A procedure based on 
interpolating primitive variables at each refinement level 
was also tried. This scheme eliminates inaccurate 
pressures at new nodes within the boundary layer but 
shows minimal improvements on convergence rates. 
Interpolation techniques are being investigated to 
eliminate this inaccuracy and help increase 
convergence rates. 


CONCLUDING REMARKS 


An adaptive finite element procedure that uses 
both triangles and quadrilateral elements is coupled 
with a multi-step Galerkin-Runge Kutta algorithm to 
model laminar hypersonic compressible flows. The finite 
element mesh is adapted in regions of high gradients by 
the use of inviscid and viscous refinement indicators. 
The facility to add in layers of structured quadrilateral 
elements at the wall for better representation of the 
boundary layer is included within the framework of the 
program. 

The ability of the adaptive finite element 
procedure to capture details of compressible flow with 
strong-inviscid viscous interactions is demonstrated by 
modelling the Mach 11.7 flow over a ramp. To ensure 
adequate resolution of complex flow features, the flow 
domain is split up into three sub-regions: (a) the sharp 
leading edge, (b) the flat plate region, and (c) the ramp 
section. 

The accuracy of the finite element solution at the 
sharp leading edge section is evaluated by comparison 
with predictions of strong interaction theories. The 
complete theory of Bertram and Blackstock yield 
estimates for surface pressure distributions while the 
heat-transfer coefficient is obtained from the strong 
interaction theory of Li and Nagamatsu. Good 
correlations between the finite element predictions and 
the strong interaction theories are obtained for the 
leading edge section. 

The finite element mesh for the ramp 
section captures the physics of the complicated viscous- 
inviscid interactions at the corner. Surface coefficients of 
heat transfer, skin friction and pressure obtained from 
the finite element calculations compare well with the 
experimental measurements of Holden. The effect of the 
leading edge shock on the surface is seen by the peak 
in heating rates up near the exit section of the ramp. 

Smooth mesh transitions within the boundary 
layer is seen to be important to model the flow behavior 
without any localized ’wiggles" in the solution. 
Elimination of all explicit dissipation within the boundary 
layer is essential to obtain accurate predictions of heat 
transfer, skin friction and pressure coefficients. The 
refinement indicator used to key in on the boundary 
layer needs to be modified to ignore enrichment at 
regions where the velocity profiles are linear. The ability 
to pick out elements that lie in locations of flow reversals 
is also needed for flow domains that include separation 
and reattachment regions. A refinement procedure that 
puts in an excessive number of elements at the surface 


is to be avoided since this leads to increased 
computational effort and slower convergence rates. 

The study highlights some of the complexities 
involved in modelling compressible viscous flow 
problems, especially, high speed flows with strong 
shock-boundary layer interactions. The physics of the 
flow as well as surface coefficients obtained from the 
adaptive finite element approach compare favorably 
with theoretical and experimental data. This trend is 
encouraging and the procedure indicates good potential 
for accurate aerotherma! load predictions in the design 
of high speed vehicles. 
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(a) Triangular enrichment 



(b) Quad enrichment 



(c) Quad enrichment with triangular 
transition 


Fig. 1 Mesh Enrichment using quadrilateral and 
triangular elements. 
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(a) Original Mesh 



(b) Refinement of elements B, C, and D 



(c) Derefinement of element C 

Fig. 3 Example of adaptation procedure using 
triangular transition elements. 



Fig. 5 Subdivision of hypersonic flow domain. 



Fig. 6 Skeleton mesh for sharp leading edge. 
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(a) Finite element mesh 

Fig. 7 Finite element mesh and density contours 
for sharp leading edge. 


Fig. 4 Hypersonic flow over a 15° ramp. 
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(b) Density contours 

Fig. 7 Finite element mesh and density contours 
for sharp leading edge. 
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Fig. 9 Pressure distribution in the vicinity of the 
leading edge. 
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(a) Finite element mesh 


Fig. 10 Finite element mesh and pressure 
distribution at the exit plane of sharp 
leading edge. 
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(b) Density contours 


Fig. 8 Finite element mesh in the vicinity of the 
leading edge. 

Fig. 1 1 Pressure distribution at exit plane of the 
leading edge section. 
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Fig. 12 Density distribution at exit plane of the 
leading edge section. 


Fig. 14 Hypersonic flow detail at sharp leading 
edge. 



(a) Surface pressure ratio 



(b) Heat transfer coefficient 


Fig. 13 Comparison of pressure and heat transfer 
coefficients for the sharp leading edge. 


Ra 

it plate seel 

— 

ion 

i 

.5 

Tos ' ’ 












x ° -1 Plate x o .7 

Fig. 15 Skeleton mesh for flat plate section. 
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(b) Density contours 

Fig. 16 Finite element mesh and density 
contours for flat plate section. 



x = .11 x = .15 

(a) Finite element mesh at inflow 



(b) Density contours at inflow section 


Fig. 1 7 Detailed mesh at inflow for flat plate 
section. 



0 .5 1.0 1.5 2.0 2.5 3.0 


Density 

Fig. 18 Comparison of inlet and exit profiles of 
flat plate section. 
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Fig. 19 Density profiles at exit of the leading 

edge and inflow to the flat plate section. 
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Fig. 20 Effect of interpolation procedure on 
pressure profiles at inflow to flat plate 
section. 
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Fig. 31 Temperature profiles at loactions up the ramp. 



Fig. 34 Comparison of skin friction coefficient on ramp 
surface. 




Fig. 32 Comparison of heat transfer coefficients on the 
ramp surface. 


Fig. 35 Comparison of surface pressure coefficients on 
ramp surface. 


* Run 21 
— Finite Element 




Fig. 33 Comparison of heat transfer coefficients on the 
ramp surface. 
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Fig. 36 Pressure profiles at location A showing effects of 
refinement. 
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Fig. 37 Velocity profiles at location A on ramp. 
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Fig. 38 Density profiles on ramp showing effects of 
refinement. 
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Fig. 39 Pressure profiles on ramp showing effects of 
varying threshold values. 
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Fig. 40 Temporal history of pressure at typical new node 
within the boundary. 



